Lqg artificial pancreas control system and related method

ABSTRACT

The invention relates to a methods and systems for determining an insulin dosing recommendation. The invention employs Linear Quadratic methodology to determine the insulin dosing recommendation based on a patient&#39;s present physiological state, which is estimated by an adaptive filter methodology employing a dynamic model, which utilizes real-time measurements of blood glucose concentration.

CROSS-REFERENCES TO RELATED APPLICATIONS

This application claims benefit under 35 U.S.C. § 119(e) to U.S.Provisional Patent Application Ser. No. 60/936,613 filed on Jun. 21,2007, and to U.S. Provisional Patent Application Ser. No. 60/964,667filed on Aug. 14, 2007, which are hereby incorporated by reference intheir entireties.

BACKGROUND OF THE INVENTION

In health, blood glucose (BG) is tightly controlled by a hormonalnetwork that includes the gut, liver, pancreas, and brain, ensuringstable fasting BG levels (˜80-100 mg/dl) and transient postprandialglucose fluctuations. Diabetes is a combination of disorderscharacterized by absent or impaired insulin action, resulting inhyperglycemia. Intensive insulin and oral medication treatment tomaintain nearly normal levels of glycemia markedly reduces chroniccomplications in both Type 1 and Type 2 diabetes, but may cause a riskof potentially life-threatening severe hypoglycemia (SH). This SHresults from imperfect insulin replacement, which may reduce warningsymptoms and hormonal defenses. Consequently, hypoglycemia has beenidentified as the primary barrier to optimal diabetes management.

DESCRIPTION OF RELATED ART

Glucose control has been studied for more than 3 decades and widelydifferent solutions have been proposed, though it is only very recentlythat technology and algorithm have come together to enable glucosecontrol elsewhere than in the ICU of a hospital.

Information relevant to attempts to provide glucose control based onintravenous (IV) glucose measure and both positive (glucose) andnegative (insulin) control actuation can be found in the followingreferences, which are not admitted to be prior art with respect to thepresent invention by inclusion in this section:

-   -   (1) Pfeiffer E F, Thum Ch, and Clemens A H: The artificial beta        cell—A continuous control of blood sugar by external regulation        of insulin infusion (glucose controlled insulin infusion        system). Horm Metab Res 487: 339-342, 1974; and    -   (2) Clemens A H: Feedback control dynamics for glucose        controlled insulin infusion system. MedProg Technol 6: 91-98,        1979.        Both of these regulators were based on a        proportional-integral-derivative strategy (PID): the injected        insulin is proportional to the difference between a fixed plasma        glucose target and the measured plasma glucose as well as to the        rate of change of plasma glucose. However, each one of these        reference suffers from the following disadvantage, which are not        admitted to have been known in the art by inclusion in this        section:    -   (1) PID control is not “model-based.” In other words, prediction        plays no role in computing control actions; and    -   (2) PID control actions are computed merely as a simple linear        function of the error signal.

Information relevant to attempts to provide glucose control based onprediction of glucose level can be found in the following references,which are not admitted to be prior art with respect to the presentinvention by inclusion in this section:

-   -   (1) Albisser A M, Leibel B S, Ewart T G, Davidovac Z, Botz C K,        and Zingg W: An artificial endocrine pancreas. Diabetes, 23:        389-396, 1974;    -   (2) Marliss E B, Murray F T, Stokes E F, Zinman B, Nakhooda A F,        Denoga A, Leibel B S, and Albisser A M: Normalization of        glycemia in diabetics during meals with insulin and glucagon        delivery by the artificial pancreas. Diabetes 26: 663-672, 1977;    -   (3) Kraegen E W, Campbell L V, Chia Y O, Meler H, and Lazarus L:        Control of blood glucose in diabetics using an artificial        pancreas. Aust N Z J Med 7: 280-286, 1977;    -   (4) Fischer U, Jutzi E, Freyse E-J, and Salzsieder E: Derivation        and experimental proof of a new algorithm for the artificial        r-cell based on the individual analysis of the physiological        insulin-glucose relationship. Endokrinologie 71:65-75, 1978; and    -   (5) Broekhuyse H M, Nelson J D, Zinman B, and Albisser A M:        Comparison of algorithms for the closed-loop control of blood        glucose using the artificial beta cell. IEEE Trans Biomed Eng        BME-28: 678-687, 1981.

Providing glucose control based on prediction of glucose level tends tocounteract the inherent inertia of exogenous insulin compared to theendogenous hormones. However, each one of these reference suffers fromone or more of the following disadvantages, which are not admitted tohave been known in the art by inclusion in this section:

-   -   (1) Since a new control profile is computed for each CGM reading        based on predictions, Model Predictive Control (MPC) is much        more computationally demanding than PID;    -   (2) The frequency of control updates may be limited by the        computational capacity of the platform on which the control is        implemented; and    -   (3) Inaccuracy in the predictive model results in poor glucose        control.

Information relevant to attempts to provide glucose control, spanning abroader range of control theory, can be found in the followingreferences, which are not admitted to be prior art with respect to thepresent invention by inclusion in this section:

-   -   (1) Salzsieder E, Albrecht G, Fischer U, and Freyse E-J: Kinetic        modeling of the glucoregulatory system to improve insulin        therapy. IEEE Trans Biomed Eng 32: 846-855, 1985;    -   (2) Fischer U, Schenk W, Salzsieder E, Albrecht G, Abel P, and        Freyse E-J: Does physiological blood glucose control require an        adaptive strategy? IEEE Trans Biomed Eng BME-34:575-582, 1987;    -   (3) Sorensen J T: A Physiologic Model of Glucose Metabolism in        Man and its Use to Design and Assess Improved Insulin Therapies        for Diabetes, Ph.D. dissertation, Department of Chemical        Engineering, MIT, 1985;    -   (4) Swan G W: An optimal control model of diabetes mellitus.        Bull Math Bio 44: 793-808, 1982;    -   (5) Fisher M E and Teo K L: Optimal insulin infusion resulting        from a mathematical model of blood glucose dynamics. IEEE Trans        Biomed Eng 36: 479-486, 1989;    -   (6) Ollerton R L: Application of optimal control theory to        diabetes mellitus. Int J Control 50: 2503-2522, 1989;    -   (7) Fisher M E: A semi closed-loop algorithm for the control of        blood glucose levels in diabetics. IEEE Trans Biomed Eng 38:        57-61, 1991; and    -   (8) Kienitz K H and Yoneyama T: A robust controller for insulin        pumps based on H-infinity theory. IEEE Trans Biomed Eng 40:        1133-1137, 1993.        However, each one of these reference suffers from one or more of        the following disadvantages, which are not admitted to have been        known in the art by inclusion in this section:    -   (1) all were concerned with IV sensing, and IV action; and    -   (2) most relied on some approximate modeling of human        physiology.

Attempts have also been made to provide glucose control based onself-monitoring of blood glucose (SMBC) to adjust the dosing of insulindelivered via injections or insulin pump. However, such attempts toprovide glucose control suffer from many disadvantages, which are notadmitted to have been known in the art by inclusion in this section.Glucose is measured at infrequent (<5/day) and irregular times duringthe day and insulin is injected subcutaneously according to both thesemeasures and the estimated amount of carbohydrates ingested. Dependingon the treatment strategy, the insulin is either injected continuously(basal rate) or discretely (boluses) via a pump, or only discretely, viainjections containing both fast acting and long acting insulin. In bothcases therelation between the amount of insulin injected and themeasured plasma glucose is determined by the care practitioner and thepatient based on past experience and initial rule of thumbs (1800-ruleand 450-rule). Insulin boluses are traditionally calculated in twophases: First, the amount of insulin is computed that is needed by aperson to compensate for the carbohydrate content of an incoming meal.This is done by estimating the amount of carbohydrates to be ingestedand multiplying by each person's insulin/carbohydrate ratio. Second, thedistance between actual blood glucose (BG) concentration and individualtarget level is calculated and the amount of insulin to reach target thetarget is computed. This is done by multiplying the (BG—target)difference by an individual insulin correction factor. It is thereforeevident that a good assessment of each person's carbohydrate ratio andcorrection factor is critical for the optimal control of diabetes.

Attempts have been made to develop regulation systems (e.g. artificialpancreas) to control insulin delivery in people with diabetes based onobserving and acting upon the glucose/insulin levels using real-timemeasurements at sampling frequencies less than or equal to 5 minutes.Some control efforts have focused on the Subcutaneous-Subcutaneous(SC-SC) route as it is the most likely to be easily mass marketed and itrelies on readily available technologies. However, each one of thesetechnologies suffer from one or more of the following disadvantages,which are not admitted to have been known in the art by inclusion inthis section:

-   -   (1) The continuous sensors currently available experience delays        estimated between 10 and 20 minutes.    -   (2) The continuous sensors' accuracy is still lower than for        example finger stick measurement (SMBG) and therefore none of        the currently available sensors have been approved for        ‘replacement’ by the Food & Drug Administration (FDA), therefore        precluding their use as such in clinical decision.    -   (3) Subcutaneous injection of insulin imposes an additional        actuation delay, the exogenous insulin being first transported        from the injection site to the central vascular system and only        then following the pathway of exogenous IV injected insulin.

Information relevant to attempts to provide glucose control based onimplantable sensors and insulin pumps can be found in the followingreferences, which are not admitted to be prior art with respect to thepresent invention by inclusion in this section:

-   -   (1) Leblanc H, Chauvet D, Lombrail P, Robert J J: Glycemic        control with closed-loop intraperitoneal insulin in type I        diabetes. Diabetes Care, 9: 124-128, 1986;    -   (2) J L Selam, P Micossi, F L Dunn, and DM Nathan: Clinical        trial of programmable implantable insulin pump for type I        diabetes, Diabetes Care 15: 877-885, 1992;    -   (3) E Renard: Implantable closed-loop glucose-sensing and        insulin delivery: the future for insulin pump therapy, Current        Opinion in Pharmacology, 2(6), 708-716, 2002; and    -   (4) R. Hovorka: Continuous glucose monitoring and closed-loop        systems Diabetic Medicine 23 (1), 1-12, 2006.        The sensors are implantable directly into an artery, and were        believed to be closer to IV sensing, and therefore less inclined        to delays and errors. However, recent studies have shown that        these sensors suffer from delays equivalent to SC sensors, even        though the implantable sensors sample blood directly.        Additionally, these attempts suffer from one or more of the        following disadvantages, which are not admitted to have been        known in the art by inclusion in this section:    -   (1) Surgery is required to insert the technologies;    -   (2) The implanted devices have a limited life time, from 3 to 18        months; and    -   (3) Despite expectations that implantable pumps would be more        efficient than SC pumps, because they more closely mimic natural        peritoneal injection of insulin, the technologies suffer from        insulin aggregation.

Information relevant to attempts to regulate glucose homeostasis havecan be found in the following references, which are not admitted to beprior art with respect to the present invention by inclusion in thissection:

-   -   (1) El-Khatib F E, Jiang J, Damiano E R: Adaptive closed-loop        control provides blood-glucose regulation using dual        subcutaneous insulin and glucagon infusion in diabetic swine.        Diabetes Sci Tech, 1:181-192, 2007;    -   (2) Hovorka R., Canonico V., Chassin L. J., Haueter U.,        Massi-Benedetti M., Federici M. O., Pieber T. R., Schaller H.        C., Schaupp L., Vering T.: Nonlinear model predictive control of        glucose concentration in subjects with type 1 diabetes.        Phsyiological Measurement, 25 (4) 905-920, 2004; and    -   (3) Steil G M, Rebrin K, Darwin C, Hariri F, Saad M F:        Feasibility of automating insulin delivery for the treatment of        type 1 diabetes. Diabetes. December; 55(12):3344-50, 2006.

Information relevant to Linear Quadratic Gaussian (LQG) controlmethodology, can be found in the following reference, which is notadmitted to be prior art with respect to the present invention byinclusion in this section: B. D. O. Anderson and J. B. Moore. OptimalControl: Linear Quadratic Methods. Prentice-Hall, Englewood Cliffs, N.J.1989. The LQG control problem lies within a larger class of controldesign methodologies, including discrete-time Markovian control problemsand robust control. Information relevant to discrete-time Markoviancontrol problems can be found in the following references, which are notadmitted to be prior art with respect to the present invention byinclusion in this section:

-   -   (1) S. D. Patek, Partially Observed Stochastic Shortest Path        Problems with Approximate Solution by Neuro-Dynamic Programming,        IEEE Transactions on Systems, Man, and Cybernetics, in press;    -   (2) S. D. Patek, Policy Iteration Type Algorithms for Recurrent        State Markov Decision Processes, Computers and Operations        Research, 31(14), December 2004, 2333-2347, 2004;    -   (3) S. D. Patek, On Terminating Markov Decision Processes with a        Risk Averse Objective Function, Automatica, 37(9), September        2001, 1379-1386, 2001; and    -   (4) E. Campos-Nanez and S. D. Patek, Dynamically Identifying        Regenerative Cycles in Simulation-Based Optimization Algorithms        for Markov Chains, IEEE Transactions on Automatic Control,        49(4), June 2004, 1022-1025, 2004.        Information relevant to robust control problems can be found in        the following references, which are not admitted to be prior art        with respect to the present invention by inclusion in this        section:    -   (1) J. W. Kamman, S. D. Patek, and S. A. Hoeckley, Application        of Multivariable Linear Control Design to Marine Towed Systems,        AIAA Journal of Guidance Control, and Dyanmics, 19(6),        1246-1251, 1996; and    -   (2) S. D. Patek and M. Athans, Optimal Robust Ho, Control, IEEE        Conference on Decision and Control (CDC 1994), New York,        December 1994, 3539-3544, 1994.

For the foregoing reasons, there remains a need for a system to computeoptimal insulin injection dosage amounts based on continuous glucosemonitoring.

BRIEF SUMMARY OF THE INVENTION

Versions of the present invention relate to a method for determining aninsulin dosing recommendation. The method comprises applying an adaptivefilter scheme to estimate a present physiological state of a patientbased on real-time measurements of blood glucose concentration in thepatient, wherein the adaptive filter scheme models the presentphysiological state of the patient using a dynamic model. The methodfurther comprises using Linear Quadratic methodology to determine aninsulin dosing recommendation based on the present physiological stateso estimated. According to particularly preferred versions of theinvention, the insulin dosing recommendation is defined as a linearcombination of a gain vector (S) and a state vector (X), so as tominimize a quadratic cost function.

According to preferred versions of the method, the dynamic model is anaugmented minimal model of glucose kinetics, an augmented reduced mealmodel, or an augmented meal model.

According to preferred versions of the method, the dynamic model employsa plurality of physiological parameters, and the method furthercomprises an initial step of customizing the dynamic model by fixing atleast one physiological parameter.

According to other preferred versions of the method, at least onephysiological parameter is fixed at an average value, wherein theaverage value represents an average of values measured for a pluralityof subjects. In these versions, it is particularly preferred that atleast one physiological parameter is modified based on measurements ofblood glucose concentration and blood insulin concentration from thepatient to ensure absence of bias between an estimated and an actualplasma glucose concentration at steady state.

According to preferred versions of the method, the method furthercomprises determining an optimal gain for the dynamic model.

According to preferred versions of the method, the dynamic model isoptimized using a Kalman filter, or an extended Kalman filter.

According to a particularly preferred version of the present invention,the method for determining an insulin dosing recommendation comprises:taking a measurement of blood glucose concentration in a patient at eachof a plurality of actuation times; applying an adaptive filter scheme ateach actuation time to generate an estimate of the patient's presentphysiological state based on the most recent measurement of bloodglucose concentration; storing the estimate of the patient's presentphysiological state; and using Linear Quadratic methodology to determinean insulin dosing recommendation based on the estimate of the patient'spresent physiological state and based on any stored estimates of thepatient's physiological state. It is also particularly preferred thatthe actuation times are spaced at regular intervals of less than 15minutes. More preferably at regular intervals of less than 10 minutes,at regular intervals of less than 5 minutes, or at regular intervals ofless than 1 minute. Most preferably, the actuation times occurcontinuously, and the measurements of blood glucose concentration areprocessed in real-time.

The present invention also relates to an artificial pancreas controlsystem comprising an observer, and an LQ Regulator. The observerpreferably comprises a programmable computer controller, which isprogrammed to receive real-time measurements of blood glucoseconcentration in a patient, and to apply an adaptive filter scheme toestimate a present physiological state of a patient based on themeasurements of blood glucose concentration. The LQ Regulator preferablycomprises a programmable computer controller, which is programmed toreceive the present physiological state estimation from the observer, touse Linear Quadratic methodology to determine an insulin dosingrecommendation based on the present physiological state.

In preferred versions of the artificial pancreas control systemaccording to the present invention, the LQ Regulator comprises aprogrammable computer controller programmed to present the insulin doingrecommendation to the patient in an open-loop application, or directlyto an insulin pump in a closed-loop application.

BRIEF DESCRIPTION OF THE DRAWINGS

These and other features, aspects, and advantages of the presentinvention will become better understood with reference to the followingdescription and appended claims, and accompanying drawings where:

FIG. 1: shows a general schematic of the LQG based plasma glucoseregulator;

FIG. 2: shows a graph for choosing an optimal ‘fit all’ Q, i.e., findingthe optimal LQ weight for a population, based on simulation and the timeto target and time below target criteria;

FIG. 3: shows a diagram of the target tracking extension of the LQregulator;

FIG. 4: shows a hierarchical diagram of application of the LQG basedglucose regulation method, organized by level of generalization; and

FIG. 5: shows a simulation of automatic control using SC sensor with 1minute actuation time, standard Kalman Filter, predictive actuation.

DETAILED DESCRIPTION OF THE INVENTION

The present invention may be understood more readily by reference to thefollowing detailed description of preferred embodiments of the inventionas well as to the examples included therein.

All numeric values are herein assumed to be modified by the term“about,” whether or not explicitly indicated. The term “about” generallyrefers to a range of numbers that one of skill in the art would considerequivalent to the recited value (i.e., having the same function orresult). In many instances, the term “about” may include numbers thatare rounded to the nearest significant figure.

In LQG control, an actuation signal is computed to minimizesquared-error deviations from a nominal operating point, which in thiscase corresponds to tight glycemic control around a reference glucoseconcentration of 100 mg/dl. The feedback control law is derived from alinearized model of the system dynamics, which is assumed to beperturbed by Gaussian white noise disturbances and measurement errors.Then using the so-called separation principle of LQG, feedback controllaws can be computed easily as a linear combination of least-squaresestimates of the states of the system through feedback gains designed tominimize the integral of weighted squared errors and control signals ina noise-free system (i.e. LQR feedback gains).

The methods and Artificial Pancreas Control Systems according to thepresent invention can be used in a wide range of applications, describedin detail below, but arise from a common general mechanism described inFIG. 1.

As shown in FIG. 1, observer (1) receives information from sensor (2)and/or from injection device (3). This information can include, forexample, a measurement of blood glucose concentration in a patient.Observer (1) then models the present physiological state of the patient.The modeled physiological state depends on the observer type (4), theobserver parameters (5), the model structure (6), and the subjectcharacteristics (7). The estimated physiological state is then appliedto a Linear Quadratic (LQ) regulator (8). The LQ regulator (8)determines an insulin dosing recommendation, which is preferably used todetermine an appropriate insulin injection (9). The determination of theinsulin dosing recommendation depends on the actuation frequency (10),the target type (11), the model structure (12), the subjectcharacteristics (13), and the discretization mitigation (14).

According to preferred versions of the present invention, blood glucoseis measured frequently (<15 min), and these measures and informationabout insulin treatment are provided to a state observer. The stateobserver estimates the physiological state of the subject (differentmodels impose different ‘physiological states’, including plasma glucoseand insulin concentrations), while the Linear Quadratic regulator usespresent and past estimations of the state to compute the needed insulintreatment.

While some characteristics of the methods and Artificial PancreasControl Systems of the present invention are highly application specific(e.g. target type, or injection device feedback) the bulk of the methodremains the same across application fields. In this section we presentthese general components.

Preferably at each received glucose measure (from the continuous sensor)the system updates its estimation of the physiological state of thesubject, i.e. it's estimation of glucose concentration (not necessarilyequal to the measure considering sensor errors) and insulinconcentration (estimated from known injection and glucose measures). Toperform this task, it is preferable to apply adaptive filtermethodology, representing the subject state by a dynamic system(equation 1.1), where X is the physiological state and Y the outcomemeasure: in this case blood glucose concentration.

$\begin{matrix}{\quad\left\{ \begin{matrix}{{\overset{\cdot}{X}(t)} = {{A_{t}{X(t)}} + {B_{t}{u(t)}}}} \\{{Y(t)} = {{C_{t}{X(t)}} + {D_{t}{u(t)}}}}\end{matrix}\; \right.} & (1.1)\end{matrix}$

From this formalism we derive the adaptive filter scheme given byequation 1.2 where {circumflex over (X)} is the estimated state of thesystem.

{circumflex over ({dot over (X)})}(t)=A _(t) {circumflex over (X)}(t)+B_(t) u(t)−L _(t)(Y(t)−C _(t) {circumflex over (X)}(t)−D _(t)u(t))  (1.2)

The methodology for selection of the state vector (X), as well as thedynamic equation parameters (A, B, C D), and the innovation gain (L) ispresented below.

Preferably, the choice of the state vector is made via two antagonisticcriteria: the higher the dimension of the model (i.e. the number ofmodeled quantities: e.g. plasma glucose, plasma insulin, remote insulin)the more precisely the model can describe observed dynamics, and as suchestimate the modeled quantities; but high dimensions also render theestimation procedure sensitive to noise in the observed signal (BG),lowering the precision of the state estimate, and can preclude the useof subject specific regulation (see below), due to unobservable statesand non identifiable parameters. Preferred versions of the presentinvention, therefore utilize one of three dynamic models, withincreasing dimensionality. All three models are original, though attheir core lie published differential equations. The models, presentedbelow, require increasingly difficult parameters estimation; however,selection of a dynamic model preferably depends on the need of the user.

The first model is augmented Minimal Model of Glucose Kinetics (AMMGK).Presented for the first time in 1979 (See, RN Bergman, YZ Ider, CRBowden, and C Cobelli. Quantitative estimation of insulin sensitivity.Am J Physio, 236:E667-EE677, 1979) the Minimal Model of Glucose Kinetics(MMGK) has since become the most used and well-known model of glucosehomeostasis in humans. Its simplicity (the original model has 2 states)has allowed the introduction of key concepts like insulin sensitivityand remote insulin.

Preferred versions of the present invention employ augmented MMGK totake into account the transport of insulin from the injection site tothe central circulation, and the transport of glucose from the stomachto the central circulation. The model is presented in equation (1.3),differential equations 1 and 3 are directly taken from the common formof the MMGK, and insulin kinetics are added in line 4. Since preferredversions of the present invention measure glucose subcutaneously, it ispreferable to add a simple diffusion equation (line 2) to represent thetransport of glucose from the compartment of interest (plasma) to themeasuring site (interstitium). Lines 5 and 6 model the subcutaneousinsulin transport after injection (J), while line 7 models the glucoserate of appearance from a meal (D). This model leads to a state vectorof dimension 7, as shown in equation (1.3).

$\quad\begin{matrix}\left\{ \begin{matrix}{\overset{\cdot}{G} = {{- {S_{g}\left( {G - G_{b}} \right)}} - {XG} + \frac{{Ra}(t)}{V_{g}}}} \\{{\overset{\cdot}{G}}_{i} = {{- \frac{1}{\tau}}\left( {G_{i} - G} \right)}} \\{\overset{\cdot}{X} = {{{- p_{2}}X} + {p_{3}\left( {I_{p} - I_{P_{basal}}} \right)}}} \\{\overset{\cdot}{I} = {{- {nI}} + \frac{{k_{a\; 1}I_{{SQ}\; 1}} + {k_{a\; 2}I_{{SQ}\; 2}}}{V_{i}}}} \\{{\overset{\cdot}{I}}_{{SQ}\; 1} = {{{- \left( {k_{a\; 1} + k_{d}} \right)}I_{{SQ}\; 1}} + {J(t)}}} \\{{\overset{\cdot}{I}}_{{SQ}\; 2} = {{{- k_{a\; 2}}I_{{SQ}\; 2}} + {k_{d}I_{{SQ}\; 2}}}} \\{{\overset{\cdot}{R}a} = {{- \frac{1}{\tau_{meal}}}\left( {{Ra} - {D(t)}} \right)}}\end{matrix} \right. & (1.3)\end{matrix}$

The second model is the reduced Meal Model (RMM). (Dalla Man et al.presented in 2007 a much more complicated model to represent glucose andinsulin fluxes during a meal (C. Dalla Man, R. A. Rizza, and C. Cobelli.Meal simulation model of the glucose-insulin system. IEEE Trans BiomedEng, in press 2007.) Using triple tracers, each flux was separatelyidentified in more than 200 subjects and modeled via the simplestavailable equations.

Preferred versions of the present invention make use of this model byaugmenting it in much the same way as the MMGK, using two subcutaneousinsulin equations, an interstitial glucose equation and modeling themeal via a single equation. This departs from the initial model, whichincluded a strongly non linear modeling of the rate of appearance ofglucose, see next paragraph. The final result is a dynamic system ofdimension 11, as shown in equation (1.4.)

(1.4) $\left\{ \begin{matrix}{{\overset{\cdot}{G}}_{p} = {{{- \left( {k_{2} + {kp}_{2}} \right)}G_{p}} + {k_{1}G_{t}} - U_{ii} - {{kp}_{3}I_{d}} + {{Ra}(t)} + {kp}_{1}}} \\{{\overset{\cdot}{G}}_{t} = {{{- k_{1}}G_{t}} + {k_{2}G_{p}} - {\frac{{Vm}_{0} + {{Vm}_{X}X}}{{Km}_{0} + G_{t}}G_{t}}}} \\{{\overset{\cdot}{G}}_{i} = {{- \frac{1}{\tau_{IG}}}\left( {G_{i} - \frac{G_{p}}{V_{g}}} \right)}} \\{{\overset{\cdot}{I}}_{d} = {- {k_{i}\left( {I_{d} - I_{1}} \right)}}} \\{{\overset{\cdot}{I}}_{1} = {- {k_{i}\left( {I_{1} - \frac{I_{p}}{V_{i}}} \right)}}} \\{{\overset{\cdot}{I}}_{p} = {{{- \left( {m_{2} + m_{4}} \right)}I_{p}} + {m_{1}I_{l}} + {k_{a\; 1}I_{{SQ}\; 1}} + {k_{a\; 2}I_{{SQ}\; 2}}}} \\{{\overset{\cdot}{I}}_{l} = {{{- \left( {m_{1} + m_{3}} \right)}I_{p}} + {m_{2}I_{p}}}} \\{\overset{\cdot}{X} = {- {p_{2u}\left( {X - \frac{I_{p}}{V_{i}} + I_{b}} \right)}}} \\{{\overset{\cdot}{I}}_{{SQ}\; 1} = {{{- \left( {k_{a\; 1} + k_{d}} \right)}I_{{SQ}\; 1}} + {J(t)}}} \\{{\overset{\cdot}{I}}_{{SQ}\; 2} = {{{- k_{a\; 2}}I_{{SQ}\; 2}} + {k_{d}I_{{SQ}\; 2}}}} \\{{\overset{\cdot}{R}a} = {{- \frac{1}{\tau_{meal}}}\left( {{ra} - {D(t)}} \right)}}\end{matrix} \right.$

The third model is the Augmented Meal Model (AMM). The model presentedbelow is an extension of the model presented by Dalla Man, like theoriginal it includes the non linear gastro-intestinal model (3 states),the 2 glucose compartments and the 5 insulin related states. However,according to preferred versions of the present invention the model isaugmented by 2 additional insulin states for transport and an additionalglucose compartment (interstitium). The resulting state, utilized inpreferred versions of the present invention, is of dimension 13, asshown in equation (1.5).

                                  (1.5) $\left\{ {{\begin{matrix}{{\overset{\cdot}{G}}_{p} = {{{- \left( {k_{2} + {kp}_{2}} \right)}G_{p}} + {k_{1}G_{t}} - U_{ii} - {E\left( G_{p} \right)} - {{kp}_{3}I_{d}} + \frac{{fk}_{{ab}\; s}Q_{gat}}{BW} + {kp}_{1}}} \\{{\overset{\cdot}{G}}_{t} = {{{- k_{1}}G_{t}} + {k_{2}G_{p}} - {\frac{{Vm}_{0} + {{Vm}_{X}X}}{{Km}_{0} + G_{t}}G_{t}}}} \\{{\overset{\cdot}{G}}_{i} = {{- \frac{1}{\tau_{IG}}}\left( {G_{i} - \frac{G_{p}}{V_{g}}} \right)}} \\{{\overset{\cdot}{I}}_{d} = {- {k_{i}\left( {I_{d} - I_{1}} \right)}}} \\{{\overset{\cdot}{I}}_{1} = {- {k_{i}\left( {I_{1} - \frac{I_{p}}{V_{i}}} \right)}}} \\{{\overset{\cdot}{I}}_{p} = {{{- \left( {m_{2} + m_{4}} \right)}I_{p}} + {m_{1}I_{l}} + {k_{a\; 1}I_{{SQ}\; 1}} + {k_{a\; 2}I_{{SQ}\; 2}}}} \\{{\overset{\cdot}{I}}_{l} = {{{- \left( {m_{1} + m_{3}} \right)}I_{p}} + {m_{2}I_{p}}}} \\{\overset{\cdot}{X} = {- {p_{2u}\left( {X - \frac{I_{p}}{V_{i}} + I_{b}} \right)}}} \\{{\overset{\cdot}{I}}_{{SQ}\; 1} = {{{- \left( {k_{a\; 1} + k_{d}} \right)}I_{{SQ}\; 1}} + {J(t)}}} \\{{\overset{\cdot}{I}}_{{SQ}\; 2} = {{{- k_{a\; 2}}I_{{SQ}\; 2}} + {k_{d}I_{{SQ}\; 2}}}} \\{{\overset{\cdot}{Q}}_{{sto}\; 1} = {{{- k_{gri}}Q_{{sto}\; 1}} + {D(t)}}} \\{{\overset{\cdot}{Q}}_{{sto}\; 2} = {{{- {k_{empt}\left( {Q_{{sto}\; 1} + Q_{{sto}\; 2}} \right)}}Q_{{sto}\; 2}} + {k_{gri}Q_{{sto}\; 1}}}} \\{{\overset{\cdot}{Q}}_{gut} = {{{- k_{{ab}\; s}}Q_{gut}} + {{k_{empt}\left( {Q_{{sto}\; 1} + Q_{{sto}\; 2}} \right)}Q_{{sto}\; 2}}}}\end{matrix}\mspace{79mu}{Where}\mspace{14mu} E\left( G_{p} \right)} = \left\{ {{\begin{matrix}{{ke}_{1}\left( {G_{p} - {ke}_{2}} \right)} & {{{if}\mspace{14mu} G_{p}} > {ke}_{2}} \\0 & {otherwise}\end{matrix}\mspace{79mu}{and}\mspace{14mu}{k_{empt}\left( {Q_{{sto}\; 1} + Q_{{sro}\; 2}} \right)}} = {{k_{m\; i\; n} + {\frac{k_{m\;{ax}} - k_{m\; i\; n}}{2}\left( {2 + {\tanh\left( {{aa}\left( {Q_{{sto}\; 1} + Q_{{sto}\; 2} - {b\mspace{14mu}{dose}}} \right)} \right)} + {\tanh\left( {{cc}\left( {Q_{{sto}\; 1} + Q_{{sto}\; 2} - {d\mspace{14mu}{dose}}} \right)} \right)}} \right)\mspace{79mu}{aa}}} = {{\frac{2.5\mspace{14mu}{dose}}{1 - b}\mspace{79mu}{cc}} = \frac{2.5\mspace{14mu}{dose}}{d}}}} \right.} \right.$

The state estimation filter presented in equation (1.2) relies onparameters of different nature to be usable: first the matrices A_(t),B_(t), C_(t), and D_(t) are derived using the chosen model (see above),the type of fit to the subject, the insulin used and the optimizationprocedure (see below). Then, depending on the characteristics of theprocess noise and the sensor noise, the optimization technique is usedto compute the filter gain L_(t). In the next paragraphs we describe indetail how to select preferable parameters in different scenarios.

As subjects can be widely different from one another, it is preferable,in order to observe a state vector, to tailor the used model to aparticular subject. As preferred versions of the present inventionemploy ‘physiological’ models, some parameters are readily available,e.g. body weight (BW). However, most parameters are highly modelspecific and require some data collection before the regulator can beput online. In the following section we describe three preferred choicesfor the model parameters estimation, ranging from a fixed parameters,plug 'n play type, methodology to a fully subject specific observer.

According to preferred versions of the present invention, parameters canbe fixed at population values. Each model has been fitted to enoughsubjects to derive averages of the parameters. Using the averages of thepopulation, it is preferable to construct a model that is subjectindependent, and therefore does not require any tuning before use.Particularly preferred values of the parameters are reported below:

AMMGK: Sg=0.0094; Vg=2.5; p2=0.0265; p3=6.724e-6; Gb=142; n=0.2;Vi=0.125; m2=0.3616; m3=0.306; m4=0.1446; ka1=0.002; ka2=0.0211;kd=0.0166; p2u=0.0276; tauIG=0.2; Ib=104; tauMEAL=0.055;

RMM: Vg=1.834; k1=0.0702; k2=0.1151; kp2=0.0061; Uii=1; kp3=0.0087;kp1=5.1207; Vm0=5.3263; Vmx=0.0417; Km0=234.0043; ki=0.0075; Vi=0.0503;m2=0.3616; m3=0.306; m4=0.1446; ka1=0.002; ka2=0.0211; kd=0.0166;p2u=0.0276; tauIG=0.2; Ib=104;

AMM: Vg=1.834; k1=0.0702; k2=0.1151; kp2=0.0061; Uii=1; kp3=0.0087;kp1=5.1207; Vm0=5.3263; Vmx=0.0417; Km0=234.0043; ki=0.0075; Vi=0.0503;m2=0.3616; m3=0.306; m4=0.1446; ka1=0.002; ka2=0.0211; kd=0.0166;p2u=0.0276; kabs=0.2386; kmax=0.0384; kmin=0.0089; b=0.7705; f=0.9;d=0.1714; ke1=0.0005; ke2=339; tauIG=0.2; Ib=104; dose=90;

Mitigated population parameters for pump users: While populationparameters are easy to use they suffer from a major drawback: at rest(steady state) there exists a bias between the estimated and the actualplasma glucose concentration, therefore making any regulator biased aswell. To ensure the absence of bias while avoiding the full model fit itis particularly preferable to employ a methodology based on a singleblood draw:

For AMM and RMM, a preferred procedure is as follows: (1) Draw blood(quantity depending on the assay used) after a minimum of 4 hoursfasting and without insulin bolus or changes in insulin basal rate.Measure blood glucose concentration [mg/dl] and blood insulinconcentration [pmol/L]. Record basal rate [mU/min]. (2) Modify thepopulation parameters as follows: Vi=2.7648*basal_rate/Iss;Vg=349.5685/Gss. This modification cancels the bias by adjusting thediffusion volumes so that the population steady state quantities, whichresult from the population parameters and therefore are constant in thiscase (Ipss=2.7648*basal_rate pmol·kg⁻¹; Gpss=349.5685 mg·kg⁻¹)correspond to the measured concentration (Iss and Gss).

For AMMGK, it is preferable to follow the same procedure as before,replacing Ib and Gb respectively by Iss and Gss, and replacing Vi by5*basal_rate/Iss.

To achieve a fully subject specific model particularly preferredversions of the present invention estimate most of the model parameters.Such estimations allow the estimator not only to be bias free at steadystate, like the previous method, but to also capture all the dynamics ofthe system, allowing for a bias free estimation both in transient mode(when glucose is not stabilized at a constant value, which is most ofthe time) and steady state. This estimation is relatively easy for theAMMGK, more cumbersome for the RMM, and extremely difficult for the fullMM.

Preferred versions of the present invention are based on analysis of anOGTT or IVGTT, which can be performed in a doctor's office. (1) Thesubject presents himself/herself after fasting for 8 hours having notchanged the insulin basal rate for 2 hours. (2) 0.5 g/kg⁻¹ of bodyweight of glucose is ingested (OGTT) or injected intravenously (IVGTT).(3) Inject (subcutaneously)×unit of insulin, x being calculated as perthe patient carbohydrate ratio. (4) Venous blood is drawn and plasmaglucose and insulin concentration are measured at minutes 0 5 10 30 4560 90 120. (5) Using the driving function methodology, and the weightedleast square criterion adjust the model parameters: the parameters arefound to minimize the weighted sum of square of the difference betweenmodel and measure at each sampling point. Some parameters are estimatedusing the insulin measures, others the glucose. It is not advised toestimate all the parameters in one run.

Insulin characteristics: In all three models the transport of insulinfrom interstitium to central circulation is preferably modeled via 2differential equations depending on 3 parameters. These 3 parametersdepend of course on the subject but can also change with different typeof insulin or mixtures thereof.

Versions of the present invention utilize fixed parameters. In theseversions, all insulin are considered have the same dynamics, coarsely,and ka1, ka2 and kd are set to population values, or subject specificvalues depending on method chose.

Preferred versions of the present invention employ insulin specificparameters. In these versions, based on response curve provided by theindustry, each insulin types is fit with a specific set of parameters,using the driving function methodology described above.

Preferred versions of the present invention employ optimal innovationgain selection. It is to be noted that while the proposed models are nonlinear, the estimator and the dynamic model it is based on are linearstate space equations. Therefore, preferred versions of the presentinvention not only to determine an optimal gain (L), but also determinethe matrices A_(t), B_(t), C_(t), and D_(t). In the next paragraphs wepresent two preferred optimization methods: the Kalman filter and theextended kalman filter.

Kalman filter optimization: In the Kalman methodology the state equationmatrices (A, B, C, and D) are constant and correspond to thelinearization of the physiological model around a specified operationpoint. The operation point should be chosen as the steady state targetof the regulator, or if not available (control between bounds) as theexpected average of the state.

Linearization is done as follows:

$\begin{matrix}{\mspace{76mu}{{AMMGK}:}} & \; \\{\mspace{79mu}{\zeta_{op} = {{\begin{bmatrix}G_{target} \\G_{target} \\{S_{g}\frac{G_{target} - G_{b}}{G_{target}}} \\{{S_{g}\frac{{p_{2}\left( {G_{target} - G} \right)}_{b}}{p_{3}G_{target}}} + I_{b}} \\\frac{\left( {\frac{S_{g}{p_{2}\left( {G_{target} - G} \right)}_{b}}{p_{3}G_{target}} + I_{b}} \right)}{\left( {k_{a\; 1} + k_{d}} \right){nV}_{i}} \\\frac{k_{d}\left( {\frac{S_{g}{p_{2}\left( {G_{target} - G} \right)}_{b}}{p_{3}G_{target}} + I_{b}} \right)}{{k_{a\; 2}\left( {k_{a\; 1} + k_{d}} \right)}{nV}_{i}} \\0\end{bmatrix}\mspace{20mu} A} = \begin{bmatrix}{{- S_{g}} - X_{op}} & 0 & G_{op} & 0 & 0 & 0 & \frac{1}{V_{g}} \\\frac{1}{\tau} & \frac{1}{\tau} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & {- p_{2}} & p_{3} & 0 & 0 & 0 \\0 & 0 & 0 & {- n} & \frac{k_{a\; 1}}{V_{i}} & \frac{k_{a\; 2}}{V_{i}} & 0 \\0 & 0 & 0 & 0 & {- \left( {k_{a\; 1} + k_{d}} \right)} & 0 & 0 \\0 & 0 & 0 & 0 & k_{d} & {- k_{a\; 2}} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & {- \frac{1}{\tau_{meal}}}\end{bmatrix}}}\;} & (1.6) \\{\mspace{79mu}{B = {{\begin{bmatrix}0 \\0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}\mspace{20mu} G} = {{\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\\frac{1}{\tau_{meal}}\end{bmatrix}\mspace{20mu} C} = {{\begin{bmatrix}0 \\1 \\0 \\0 \\0 \\0 \\0\end{bmatrix}\mspace{25mu} D} = \begin{bmatrix}0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0\end{bmatrix}}}}}\;} & \; \\{\mspace{79mu}{{RMM}:}} & \; \\{\mspace{79mu}{\zeta_{op} = {{\begin{bmatrix}{G_{pop} = {V_{g}G_{target}}} \\{G_{top} = \frac{{- \beta} \pm \sqrt{\beta^{2} - {4{\alpha\gamma}}}}{2\alpha}} \\G_{target} \\{I_{dop} = \frac{{{- \left( {k_{1} + k_{2}} \right)}G_{pop}} + {k_{2}G_{top}} + {{kp}\; 1} - {Uii}}{{kp}_{3}}} \\{I_{1\;{op}} = I_{dop}} \\{I_{pop} - \frac{I_{dop}}{V_{i}}} \\{I_{lop} = \frac{m_{2}I_{pop}}{m_{1} + m_{3}}} \\{X_{op} = {I_{dop} - I_{b}}} \\{I_{{sq}\; 1{op}} = \frac{\left( {m_{2} + m_{4}} \right)I_{pop}\mspace{14mu} m_{1}I_{lop}}{k_{a\; 1} + k_{d}}} \\{I_{{sq}\; 2{op}} - \frac{k_{d}I_{{sq}\; 1{op}}}{k_{a\; 2}}} \\0\end{bmatrix}\mspace{14mu}\alpha} = {{- k_{2}} - \frac{{Vm}_{X}*k_{2}}{{kp}_{3}}}}}} & (1.7) \\{\mspace{76mu}{{{where}\mspace{14mu}\beta} - {k_{1}G_{pop}} - {k_{2}{Km}_{0}} - {Vm}_{0} + {{Vm}_{X}I_{b}} - \frac{{{Vm}_{X}\left( {k_{1} + {kp}_{2}} \right)G_{pop}} - {{Vm}_{X}{kp}_{1}} + {{Vm}_{X}U_{ii}}}{{kp}\; 3}}} & \; \\{\mspace{70mu}{\gamma = {k_{1}G_{pop}{Km}_{0}}}} & \; \\{\mspace{70mu}{B - {\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}\mspace{14mu} G} - \begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\\frac{1}{\tau_{meal}}\end{bmatrix}}} & \; \\{A = \begin{bmatrix}{- \left( {k_{1} + {kp}_{2}} \right)} & k_{2} & 0 & {- {kp}_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\k_{1} & \begin{matrix}{{{- k}\; 2} -} \\\frac{\left( {{Vm}_{0} - {{Vm}_{X}X_{op}}} \right){Km}_{0}}{\left( {{Km}_{0} + G_{top}} \right)^{2}}\end{matrix} & 0 & 0 & 0 & 0 & 0 & {- \frac{{Vm}_{X}G_{top}}{{Kmo} + G_{top}}} & 0 & 0 & 0 \\\frac{1}{\tau_{IG}V_{g}} & 0 & {- \frac{1}{\tau_{IG}}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & {- k_{i}} & k_{i} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & {- k_{i}} & \frac{k_{i}}{V_{i}} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & {- \left( {m_{2} + m_{4}} \right)} & m_{1} & k_{a\; 1} & k_{a\; 1} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & {m\; 2} & {- \left( {m_{1} + m_{3}} \right)} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \frac{p_{2\; u}}{V_{i}} & 0 & {- p_{2u}} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \left( {k_{a\; 1} + k_{d}} \right)} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & k_{d} & {- k_{a\; 2}} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \frac{1}{\tau_{meal}}}\end{bmatrix}} & \; \\{\mspace{79mu}{C = {{\begin{bmatrix}0 \\0 \\1 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0\end{bmatrix}\mspace{14mu} D} = \begin{bmatrix}0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0\end{bmatrix}}}} & \; \\{\mspace{79mu}{{AMM}:}} & \; \\{\mspace{79mu}{{\zeta_{op} - {\begin{bmatrix}{G_{pop} = {V_{g}G_{target}}} \\{G_{top} = \frac{{- \beta} \pm \sqrt{\beta^{2} - {4{\alpha\gamma}}}}{2\alpha}} \\G_{target} \\{I_{dop} = \frac{{{- \left( {k_{1} + k_{2}} \right)}G_{pop}} + {k_{2}G_{top}} + {{kp}\; 1} - {Uii}}{{kp}_{3}}} \\{I_{1{op}} = I_{dop}} \\{I_{pop} = \frac{I_{dop}}{V_{i}}} \\{I_{lop} = \frac{m_{2}I_{pop}}{m_{1} + m_{3}}} \\{X_{op} = {I_{dop} - I_{b}}} \\{I_{{sq}\; 1{op}} = \frac{{\left( {m_{2} + m_{4}} \right)I_{pop}} - {m_{1}I_{lop}}}{k_{a\; 1} + k_{d}}} \\{I_{{oq}\; 2{op}} = \frac{k_{\overset{\cdot}{a}}I_{{sq}\; 1\;{op}}}{k_{a\; 2}}} \\0 \\0 \\0 \\0\end{bmatrix}\mspace{14mu}\alpha}} = {{- k_{2}} - \frac{{Vm}_{X}k_{2}}{{kp}_{3}}}}} & (1.8) \\{\mspace{79mu}{{{where}\mspace{14mu}\beta} - {k_{1}G_{pop}} - {k_{2}{Km}_{0}} - {Vm}_{0} + {{Vm}_{X}I_{b}} + \frac{{{Vm}_{X}\left( {k_{1} + {kp}_{2}} \right)G_{pop}} - {{Vm}_{X}{kp}_{1}} + {{Vm}_{X}U_{ii}}}{{kp}\; 3}}} & \; \\{\mspace{79mu}{\gamma = {k_{1}G_{pop}{Km}_{0}}}} & \; \\{\mspace{76mu}{B - {\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\1 \\0 \\0 \\0 \\0\end{bmatrix}\mspace{14mu} G} - \begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}}} & \; \\{A = \begin{bmatrix}{- \left( {k_{1} + {kp}_{2}} \right)} & k_{2} & 0 & {- {kp}_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{{fk}_{{ab}\; s}}{BW} \\k_{1} & \begin{matrix}{{{- k}\; 2} -} \\\frac{\begin{matrix}\left( {{Vm}_{0} +} \right. \\{\left. {{Vm}_{X}X_{op}} \right){Km}_{0}}\end{matrix}}{\left( {{Km}_{0} + G_{top}} \right)^{2}}\end{matrix} & 0 & 0 & 0 & 0 & 0 & {- \frac{{Vm}_{X}G_{top}}{{Kmo} + G_{top}}} & 0 & 0 & 0 & 0 & 0 \\\frac{1}{\tau_{IG}V_{g}} & 0 & {- \frac{1}{\tau_{IG}}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & {- k_{i}} & k_{i} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & {- k_{i}} & \frac{k_{i}}{V_{i}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \begin{matrix}{- \left( {m_{2} -} \right.} \\\left. m_{4} \right)\end{matrix} & m_{1} & 0 & k_{a\; 1} & k_{a\; 2} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & m_{2} & \begin{matrix}{- \left( {m_{1} +} \right.} \\\left. m_{3} \right)\end{matrix} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \frac{p_{2x}}{V_{i}} & 0 & {- p_{2u}} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \begin{matrix}\left( {k_{a\; 1} +} \right. \\\left. k_{d} \right)\end{matrix} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & k_{d} & {- k_{a\; 2}} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- k_{gri}} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & k_{gri} & {- k_{emptop}} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & k_{emptop} & {- k_{{ab}\; s}}\end{bmatrix}} & \; \\{\mspace{76mu}{C = {{\begin{bmatrix}0 \\0 \\1 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0\end{bmatrix}\mspace{14mu} D} = \begin{bmatrix}0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0\end{bmatrix}}}} & \; \\{\mspace{70mu}{{{Where}\mspace{14mu} k_{empop}} = {\frac{{k\mspace{14mu}\max} - {k\mspace{14mu}\min}}{2}\left( {2 - {\tanh^{2}\left( {{aa} \cdot b \cdot {dose}} \right)} - {\tanh^{2}\left( {{cc} \cdot d \cdot {dose}} \right)}} \right.}}} & \; \\{\mspace{70mu}{{aa} = {{\frac{2.5\mspace{14mu}{dose}}{1 - b}\mspace{14mu}{cc}} = \frac{2.5\mspace{14mu}{dose}}{d}}}} & \;\end{matrix}$

Computation of {circumflex over (X)}: Using equations (1.1) and one ofthe three linear models above, it is preferable to form the systempresented in equation (1.9). The initial equation is preferably modifiedto take into account the controlled aspect of the model, namelydisturbances such as meals are left unknown (represented as W_(n)) whileinsulin injections are considered known (U_(n)). Also, noise ispreferably added to the output (V_(n)) and the system is presented inits discrete form.

ξ_(n) =Aξ _(n-1) +BU _(n-1) +GW _(n-1)

Y _(n) =Cξ _(n) +V _(n)

Where ξ=X−ζ _(op)  (1.9)

Assuming W and V (process noise/measurement noise) have means 0 andcovariance matrices K and R respectively, it is preferable to build theestimator as:

$\begin{matrix}{\mspace{79mu}{{time}\mspace{14mu}{update}\mspace{14mu}({predict})\left\{ {{\begin{matrix}{{\hat{\xi}}_{n|{n - 1}} = {{A\;{\hat{\xi}}_{n - 1}} + {BU}_{n - 1}}} \\{P_{n}^{n - 1} = {{{AP}_{n - 1}A^{T}} + Q}} \\{L_{n} = {P_{n}^{n - 1}{C^{T}\left( {{{CP}_{n}^{n - 1}C^{T}} + R} \right)}^{- 1}}}\end{matrix}{Measurement}\mspace{14mu}{update}\mspace{14mu}({correct})} = \left\{ \begin{matrix}{{\hat{\xi}}_{n} = {{\hat{\xi}}_{n|{n - 1}} + {L_{n}\left( {Y_{n} - {C\;{\hat{\xi}}_{n|{n - 1}}}} \right)}}} \\{P_{n} = {\left( {I - {L_{n}C}} \right)P_{n}^{n - 1}}}\end{matrix}\; \right.} \right.}} & (1.10)\end{matrix}$

The choices of K and R are important, and particularly preferredversions of the present invention use K=0.01, and R=1. {circumflex over(ξ)} is our state estimate and P_(n) is the estimation error covariancematrix. So, according to preferred versions of the present invention, wehave:

{circumflex over (X)}={circumflex over (ξ)}+ζ _(op) P _(n)=E[(X−{circumflex over (X)})(X−{circumflex over (X)})^(T)]  (1.9)

Extended Kalman filter optimization: In the extended Kalman methodologythe linearization is not done off-line but within the recursive filter:the EKF takes into account the non-linear nature of the system andlinearizes for each iterations around the predicted state (see Kalmanfilter above), therefore creating A, B, C, and D matrices that changewith time. Preferred versions of the present invention derive again thelinearization as for the KF but this time we consider {circumflex over(X)}_(n|n−1) known, which will be consistent with the algorithmdescribed below.

$\begin{matrix}{\mspace{76mu}{{AMMGK}:}} & \; \\{\mspace{76mu}{\begin{bmatrix}G_{op} \\G_{iop} \\X_{op} \\I_{op} \\I_{{SQ}\; 1{op}} \\I_{{SQ}\; 2\;{op}} \\{Ra}_{op}\end{bmatrix} = {{{\hat{X}}_{n|{n - 1}}\mspace{25mu} A} = \begin{bmatrix}{{- S_{g}} - X_{op}} & 0 & G_{op} & 0 & 0 & 0 & \frac{1}{V_{g}} \\\frac{1}{\tau} & {- \frac{1}{\tau}} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & {- p_{2}} & p_{3} & 0 & 0 & 0 \\0 & 0 & 0 & {- n} & \frac{k_{a\; 1}}{V_{i}} & \frac{k_{a\; 2}}{V_{i}} & 0 \\0 & 0 & 0 & 0 & {- \left( {k_{a\; 1} + k_{d}} \right)} & 0 & 0 \\0 & 0 & 0 & 0 & k_{d} & {- k_{a\; 2}} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & {- \frac{1}{\tau_{meal}}}\end{bmatrix}}}} & (1.11) \\{\mspace{76mu}{B = {{\begin{bmatrix}0 \\0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}\mspace{11mu} G} = {{\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\\frac{1}{\tau_{meal}}\end{bmatrix}\mspace{11mu} C} = {{\begin{bmatrix}0 \\1 \\0 \\0 \\0 \\0 \\0\end{bmatrix}\; D} = \begin{bmatrix}0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0\end{bmatrix}}}}}} & \; \\{\mspace{76mu}{{RMM}:}} & \; \\{\mspace{76mu}{\begin{bmatrix}G_{pop} \\G_{top} \\G_{iop} \\I_{dop} \\I_{1{op}} \\I_{pop} \\I_{lop} \\X_{op} \\I_{{sq}\; 1\;{op}} \\I_{{sq}\; 2{op}} \\{Ra}_{op}\end{bmatrix} = {{{\hat{X}}_{n|{n - 1}}\mspace{31mu} B} = {{\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}\mspace{14mu} G} = \;{{\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\\frac{1}{\tau_{meal}}\end{bmatrix}\mspace{14mu} C} = {{\begin{bmatrix}0 \\0 \\1 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0\end{bmatrix}\mspace{25mu} D} = \begin{bmatrix}0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0\end{bmatrix}}}}}}} & (1.12) \\{A = \begin{bmatrix}{- \left( {k_{1} + {kp}_{2}} \right)} & k_{2} & 0 & {- {kp}_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\k_{1} & \begin{matrix}{{{- k}\; 2} -} \\\frac{\left( {{Vm}_{0} + {{Vm}_{X}X_{op}}} \right){Km}_{0}}{\left( {{Km}_{0} + G_{iop}} \right)^{2}}\end{matrix} & 0 & 0 & 0 & 0 & 0 & {- \frac{{Vm}_{X}G_{top}}{{Kmo} + G_{top}}} & 0 & 0 & 0 \\\frac{1}{\tau_{IG}V_{g}} & 0 & {- \frac{1}{\tau_{IG}}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & {- k_{i}} & k_{i} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & {- k_{i}} & \frac{k_{i}}{V_{i}} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \begin{matrix}{- \left( {m_{2} +} \right.} \\\left. m_{4} \right)\end{matrix} & m_{1} & k_{a\; 1} & k_{a\; 2} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & {m\; 2} & {- \left( {m_{1} + m_{3}} \right)} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \frac{p_{2u}}{V_{i}} & 0 & {- p_{2u}} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \left( {k_{a\; 1} + k_{d}} \right)} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & k_{d} & {- k_{a\; 2}} & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- \frac{1}{\tau_{meal}}}\end{bmatrix}} & \; \\{\mspace{76mu}{{AMM}:}} & \; \\{\mspace{76mu}{\begin{bmatrix}G_{pop} \\G_{lop} \\G_{iop} \\I_{dop} \\I_{1{op}} \\I_{pop} \\I_{lop} \\X_{op} \\I_{{sq}\; 1{op}} \\I_{{sq}\; 2{op}} \\Q_{{sto}\; 1{op}} \\Q_{{sto}\; 2{op}} \\Q_{gutop} \\\;\end{bmatrix} = {{{\hat{X}}_{n|{n - 1}}\mspace{14mu} B} = {{\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\1 \\0 \\0 \\0 \\0\end{bmatrix}\mspace{14mu} G} = {{\begin{bmatrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\1 \\0 \\0\end{bmatrix}\mspace{20mu} C} = {{\begin{bmatrix}0 \\0 \\1 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0\end{bmatrix}\mspace{14mu} D} = \begin{bmatrix}0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0 \\0 & 0\end{bmatrix}}}}}}} & (1.13) \\{A = \begin{bmatrix}{- \left( {k_{1} + {kp}_{2}} \right)} & k_{2} & 0 & {- {kp}_{3}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{{fk}_{{ab}\; s}}{BW} \\k_{1} & \begin{matrix}{{{- k}\; 2} -} \\\frac{\begin{matrix}\left( {{Vm}_{0} +} \right. \\{\left. {{Vm}_{X}X_{op}} \right){Km}_{0}}\end{matrix}}{\left( {{Km}_{0} + G_{top}} \right)^{2}}\end{matrix} & 0 & 0 & 0 & 0 & 0 & {- \frac{{Vm}_{X}G_{top}}{{Kmo} + G_{top}}} & 0 & 0 & 0 & 0 & 0 \\\frac{1}{\tau_{IG}V_{g}} & 0 & {- \frac{1}{\tau_{IG}}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & {- k_{i}} & k_{i} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & {- k_{i}} & \frac{k_{i}}{V_{i}} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \begin{matrix}{- \left( {m_{2} +} \right.} \\\left. m_{4} \right)\end{matrix} & m_{1} & 0 & k_{a\; 1} & k_{a\; 2} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & m_{2} & \begin{matrix}{- \left( {m_{1} +} \right.} \\\left. m_{3} \right)\end{matrix} & 0 & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & \frac{p_{2u}}{V_{i}} & 0 & {- p_{2u}} & 0 & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \begin{matrix}{- \left( {k_{a\; 1} +} \right.} \\\left. k_{d} \right)\end{matrix} & 0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & k_{d} & {- k_{a\; 2}} & 0 & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & {- k_{gri}} & 0 & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & k_{gri} & \left. k_{emptop}^{(1)} \right) & 0 \\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \begin{matrix}\left( {k_{emptop} -} \right. \\\left. {Q_{{sto}\; 2{op}}k_{emptop}^{(1)}} \right)\end{matrix} & {- k_{{ab}\; s}}\end{bmatrix}} & \; \\{\mspace{76mu}{{{Where}\mspace{14mu} k_{emptop}} = {k_{m\; i\; n} + {\frac{{k\mspace{11mu}\max} - {k\mspace{11mu}\min}}{2}\left( {2 + {\tanh\left( {{aa}\left( {Q_{{sto}\; 1{op}} + Q_{{sto}\; 2{op}} - {b \cdot {dose}}} \right)} \right)} + {\tanh\left( {{cc}\left( {Q_{{sto}\; 1{op}} + Q_{{sto}\; 2{op}} - {d \cdot {dose}}} \right)} \right)}} \right)}}}} & \; \\{\mspace{76mu}{k_{emptop} = {\frac{{k\mspace{11mu}\max} - {k\mspace{11mu}\min}}{2}\left( {2 - {\tanh^{2}\left( {{aa}\left( {Q_{{sto}\; 1{op}} + Q_{{sto}\; 2{op}} - {b \cdot {dose}}} \right)} \right)} - {{cc}\mspace{11mu}{\tanh^{2}\left( {{cc}\left( {Q_{{sto}\; 1{op}} + Q_{{sto}\; 2{op}} - {d \cdot {dose}}} \right)} \right)}}} \right)}}} & \; \\{\mspace{76mu}{{aa} = {{\frac{2.5\mspace{14mu}{dose}}{1 - b}\mspace{14mu}{cc}} = \frac{2.5\mspace{14mu}{dose}}{d}}}} & \;\end{matrix}$

Computation of {circumflex over (X)}: Using equations (1.1) and one ofthe three original models, it is preferable to form the system presentedin equation (1.4).

X _(n) =f(X _(n-1) ,U _(n-1) ,W _(n-1))

Y _(n) =h(X _(n) ,V _(n))  (1.14)

Where f and h are two non-linear functions defined by the chosen model.Assuming W and V (process noise/measurement noise) have means 0 andcovariance matrices K and R respectively. It is preferable to build theestimator as:

$\mspace{79mu}{{time}\mspace{14mu}{update}\mspace{14mu}({predict})\left\{ {\begin{matrix}{{\hat{X}}_{{n\mspace{11mu} n} - 1} = {f\left( {{\hat{X}}_{n - 1},U_{n - 1},0} \right)}} \\{P_{n}^{n - 1} = {{A_{n}P_{n - 1}A_{n}^{T}} + {G_{n}Q_{k - 1}G_{n}^{T}}}}\end{matrix}{Measurement}\mspace{14mu}{update}\mspace{14mu}({correct})\left\{ {{\begin{matrix}{L_{n} = {P_{n}^{n - 1}{C_{n}^{T}\left( {{C_{n}P_{n}^{n - 1}C_{n}^{T}} + R_{n}} \right)}^{- 1}}} \\{{\hat{X}}_{n} = {{\hat{X}}_{n|{n - 1}} + {L_{n}\left( {Y_{n} - {h\left( {{\hat{X}}_{n|{n - 1}},0} \right)}} \right)}}} \\{P_{n} = {\left( {I - {L_{n}C_{n}}} \right)P_{n}^{n - 1}}}\end{matrix}\mspace{79mu}{Where}\mspace{14mu} A_{n}},G_{n},{{and}\mspace{14mu} C_{n}\mspace{14mu}{are}\mspace{14mu}{the}\mspace{14mu}{linearized}\mspace{14mu}{matrices}}} \right.} \right.}$

The choices of K and R are important, and particularly preferredversions of the present invention use K=0.01, and R=1.

LQ regulator: The second component that is common to all applicationfields is the regulator, which transforms our knowledge of thephysiological state of the system (from last section) into an insulindosing recommendation. The LQ regulator preferably performs three steps:(i) Step 1 receives a state estimation; (ii) Step 2 uses the state and apredefined gain to compute the insulin dose; (iii) Step 3 presents theinsulin dose recommendation to the patient in open-loop application, ordirectly to the insulin pump in closed-loop application.

Acquiring the estimate of the state: The estimate is preferably acquiredat each actuation time (depending of the chosen actuation frequency) byrunning the chosen observer from last actuation time, using the laststate estimates and covariance matrix as initial values, to the currentactuation time, taking into account all intermediate sensor measures.

Calculating the optimal insulin dose: To compute the insulin dose it ispreferable to use the Linear Quadratic methodology: the insulininjection is preferably defined as a linear combination of a gain vector(S) and the state vector (X), so as to minimize a quadratic costfunction.

Assuming the subject glucose and insulin dynamics can be modeled byequation (1.1), tailored using one of the linearization scheme(equations (1.6) to (1.13), preferred versions of the present inventionconstruct the following cost function to be minimized:

J=∫X ^(T) QX+U ^(T) RU  (1.15)

The LQR gain (S) minimizing J is the solution of the following Riccatiequation:

Q±A ^(T) S±SA−(SB±N)R ⁻¹(B ^(T) S±N ^(T))=0  (1.16)

Choosing Q and R: The choice of Q and R determines both the states thatare of interest (glucose vs. insulin, plasma concentration vs. othercompartment) and the aggressiveness of the regulator (how fast it willtry to achieve the selected target, therefore risking to dip under). Ina simulation analysis, it was found that only the ration of Q over R iscritical, and therefore preferred versions of the invention fix R at 1.

The choice of Q can be achieved in different ways, which are presentedin the following paragraph, but always depends on the reliability of thesensor, the frequency of actuation (therefore the value is applicationdependent), and the subject characteristics. While not all of thisinformation is always available, it is important to note that Q needs tobe tailored for every implementation of the regulator. Preferredversions of the invention use the following structure for the Q vector:

-   -   AMMGK: Q=[q q 1 1 1 1 1]    -   RMM: Q=[q q q 1 1 1 1 1 1 1]    -   AMM: Q=[q q q 1 1 1 1 1 1 1 1 1 1]

Fixed q.: In a ‘one size fits all’ approach, preferred versions of thepresent invention, use a single vector Q for all subjects. In theseversions, it is therefore important to choose a Q vector that minimizesthe risk of dip below the target while achieving reasonable control. Thefollowing methodology is preferably employed:

-   -   (1) Using a meal simulator (e.g. See, C Dalla Man, DM. Raimondo,        RA. Rizza, C Cobelli: GIM, Simulation Software of Meal        Glucose—Insulin Model, Journal of Diabetes Science and        Technology, 1(3): Page 323-330, 2007) for 100 subjects simulate        a 45 g carbohydrate meal using the chosen regulation strategy        and q between 1 and 500.    -   (2) Record the time to target (t2tgt) defined as the number of        minutes it takes each in-silico subject to come back within 10        mg/dl of the intended target.    -   (3) Record the time below target (tbtgt) defined as the number        of minutes spent below 2 mg/dl of the intended target.    -   (4) Compute J(q) as:

$\begin{matrix}{{{J(q)} = {{\frac{1}{100}{\sum\limits_{i = 1}^{100}{t\; 2\;{{tgt}(q)}_{i}}}} + {{tbtgt}(q)}_{i}}}{{{Where}\mspace{14mu} k_{empop}} = {{- {aa}}\;\frac{{k\mspace{14mu}\max} - {k\mspace{14mu}\min}}{2}\left( {2 - {\tanh^{2}\left( {{aa} \cdot b \cdot {dose}} \right)} - {\tanh^{2}\left( {{cc} \cdot d \cdot {dose}} \right)}} \right)}}{{aa} = {{\frac{2.5\mspace{14mu}{dose}}{1 - b}{cc}} = \frac{2.5\mspace{14mu}{dose}}{d}}}} & (1.17)\end{matrix}$

-   -   (5) The optimal parameter (see, FIG. 2) is defined as:

$\begin{matrix}{q^{*} = {\underset{q}{\arg\;\min}{J(q)}}} & (1.18)\end{matrix}$

Subject Specific q.: If the subject model parameters are known it ispreferable to tailor the regulator to this particular subject, thereforeimproving dramatically the performances of the regulator. To estimatethe optimal q parameter for a specific subject we propose the followingmethod is preferably employed:

-   -   (1) Using a meal simulator (e.g. C Dalla Man, DM. Raimondo, RA.        Rizza, C Cobelli: GIM, Simulation Software of Meal        Glucose—Insulin Model, Journal of Diabetes Science and        Technology, 1(3): Page 323-330, 2007) and the subject        characteristics simulate a 45 g carbohydrate meal using the        chosen regulation strategy and q equal 1, 250 and 500,        respectively named q_(min), q_(med), and q_(max).    -   (2) For each q compute tbtgt as defined above.    -   (3) If tbtgt(q_(min))>0 restart step 1 with a lower q_(min)        value.    -   (4) If tbtgt (q_(max))=0 restart step 1 with a higher q_(max)        value.    -   (5) simulate with q=q_(med)        -   If tbtgt(q)=0            -   q_(min)=q_(med)            -   q_(med)=(q_(max)−q_(min))/2        -   else if tbtgt(q)>0            -   q_(max)=q_(med)            -   q_(med)=(q_(max)−q_(min))/2    -   (6) If (q_(max)−q_(min))>1 iterate step 3    -   (7) q*=q_(min)

Target tracking method: According to preferred versions of the presentinvention, the LQ method is adapted to track an offline optimaltrajectory. In these versions, the state equations around the targetvalue at each actuation time are linearized by coupling a meal detectionor meal announcement device and a meal trajectory (computed from pastdata or stored in the device). These versions then compute S as thesolution of the new Riccati equation (since A, B, C, D and G havechanged).

FIG. 3 shows a target tracking Linear Quadratic (LQ) control systemaccording to preferred versions of the present invention. As shown inFIG. 3, a meal model reduced state observer (19) is coupled with adisturbance predictor (16), a trajectory planner (17), and an LQ tracker(20) to track offline optimal trajectory based on Continuous GlucoseMonitoring (CGM) (15). This system generates an insulin dosingrecommendation, which is preferably used to determine an appropriateinsulin dose (18).

Calculation of the injection: The LQR methodology assumes continuous orvery frequent information about the insulin injection, which is notaccessible in many implementations. To mitigate this discrepancypreferred versions of the present invention use three methods to computethe insulin injection, knowing the LQ parameters (linear model, Q, andR). For all the following methods it is assumed that the injection iscontinuous over the actuation period, in the case of non continuousactuation (like boluses) the method needs to be adapted by taking thefinal result (Inj) and multiplying by the actuation period(t_(n)−t_(n-1)) to obtain the total dose to be injected.

Instantaneous Injection: At each actuation decision time n (to bedetermined by the actuation frequency, based on application), it ispreferable to compute the product of the LQ gain (S) by the estimatedstate vector (S).

Inj=S·{circumflex over (X)} _(n)  (1.19)

Observed average injection: At each actuation decision time n, it ispreferable to compute the average theoretical injection over the lastactuation period (time since last actuation time n−1).

$\begin{matrix}{{Inj} = {\frac{1}{t_{n} - t_{n - 1}}{\int\limits_{t_{n - 1}}^{t_{n}}{{S \cdot {\hat{X}}_{t}}{dt}}}}} & (1.20)\end{matrix}$

Where {circumflex over (X)}_(t) is equal to {circumflex over (X)}_(t*),and t* is the last time of a sensor reading such that t*≤t.

Predicted average injection: At each actuation decision time n, it ispreferable to compute the average theoretical injection over the nextactuation period (time since last actuation time n−1). This means usingthe linearized system to predict the state over the injection period andcomputing the optimal injection over the period, then computing theaverage injection to use in practice, see equation (1.21).

$\begin{matrix}{{{Inj} = {\frac{1}{t_{n + 1} - t_{n}}{\int\limits_{t_{n}}^{t_{n + 1}}{{S \cdot {\overset{\sim}{X}}_{t}}{dt}}}}}{{\overset{\overset{\cdot}{\sim}}{X}}_{t} = {{A\;{\overset{\sim}{X}}_{t}} + {BU}}}{{\overset{\sim}{X}}_{t_{n}} = {\hat{X}}_{t_{n}}}} & (1.21)\end{matrix}$

Tailoring the insulin injection to the injection device: To use aregulator properly we need to take into account the characteristics ofthe hardware, i.e. maximum and minimum injection rate for pumps, maximumand minimum amount for bolus injection devices, maximum and minimumupdate frequency for pumps, minimal increment in insulin dosing for bothpump and bolus injection devices. Each limitation has influence on thetheoretical injection calculated in the above section. Preferredversions of the present invention use several methods to mitigate theseinjection hardware shortcomings:

Minimal instantaneous error: In preferred versions of the presentinvention, the LQ optimal injection is rounded to the closest availablevalue.

Inj _(final)=round[Inj/increment]×increment  (1.22)

In the case of an insulin pump, preferred versions of the presentinvention achieve better error minimization by using the duration ofinjection. Considering an actuation period of T, and possible bolusdurations of d₁< . . . <d_(p)<T we can construct the followinginjection:

$\begin{matrix}{{{Inj}_{final} = {{{round}\left\lbrack \frac{Inj}{{increment} \times d_{i}} \right\rbrack} \times {increment}\mspace{14mu}{over}\mspace{14mu} d_{i}\mspace{14mu}{minutes}}}{{{where}\mspace{14mu} i} = {\underset{{j = 1},\ldots\mspace{14mu},p}{\arg\;\min}\left( {{{{Inj} - {{round}\left\lbrack \frac{Inj}{{increment} \times d_{j}} \right\rbrack}}} \times {increment} \times d_{j}} \right)}}} & (1.23)\end{matrix}$

Minimal cumulative error: In preferred versions of the presentinvention, the new injection is calculated as the injection optimizingthe cumulative error (since beginning of operation) between LQ injectionand actual injection:

$\begin{matrix}{{{{Inj}_{final}\left( t_{n} \right)} = {{{round}\left\lbrack \frac{{Inj}_{t_{n}} + ɛ_{t_{n}}}{increment} \right\rbrack} \times {increment}}}{ɛ_{t_{n}} = {ɛ_{t_{n - 1}} + {Inj}_{t_{n}} - {{Inj}_{final}\left( t_{n} \right)}}}{ɛ_{t_{0}} = {{Inj}_{t_{0}} - {{Inj}_{final}\left( t_{0} \right)}}}} & (1.24)\end{matrix}$

Minimal cumulative error with forgetting factor: In preferred versionsof the present invention, the new injection is calculated as theinjection optimizing the cumulative weighted error since beginning ofoperation) between LQ injection and actual injection, therefore oldunresolved errors do not influence the decision.

$\begin{matrix}{{{{Inj}_{final}\left( t_{n} \right)} = {{{round}\left\lbrack \frac{{Inj}_{t_{n}} + ɛ_{t_{n}}}{increment} \right\rbrack} \times {increment}}}{ɛ_{t_{n}} = {\sum\limits_{i = 1}^{n}{e^{\gamma{({t_{i} - t_{n}})}}\left( {{Inj}_{t_{i}} - {{Inj}_{final}\left( t_{i} \right)}} \right)}}}} & (1.25)\end{matrix}$

γ is the forgetting factor, the higher γ, the less past errors influencethe new injection. Particularly preferred versions of the presentinvention use a value of γ=0.07 min⁻¹.

The methods and Artificial Pancreas Control Systems of the presentinvention can be employed in at least three different fields ofapplication of the LQ-based regulator (See, FIG. 4):

-   -   (i) computation of optimal boluses for isolated injection        insulin treatments (meal and correction boluses),    -   (ii) open loop control using an insulin pump, where the user is        advised on pre-meal boluses and basal rate changes, and    -   (iii) closed loop glucose control, where the sensor is        automatically linked with insulin pump via the algorithm, with        only safety oversight by the user.

Each application requires different settings of the method and restrictthe use of some of the modules presented above, either because theycannot be applied, or because they are sub-optimal or could be harmful.In this section we review each application field, describing possibleimplementation and module restrictions.

Bolus treatment: injection or insulin pen—type device: The maincharacteristic of this application field is the irregular infrequentactuation time. Therefore, there is a need to estimate the effect ofinsulin over a long period of time, as opposed to open or closed-loopregulation. While the state of the system can be observed almostreal-time via the previously presented observer (iterations aretriggered by frequently arriving sensor values), patients act on theirblood sugar a few times a day (commonly 3 to 5), usually at mealtime andin case of very high blood glucose. One of the characteristic of thistype of application is the combination of different versions ofsynthetic insulin in one injection (long acting and short acting). Thesecharacteristics force the following regulator characteristics: (1)observer triggered by arriving sensor data and not only by actuationtime; (2) modified injection for boluses (see calculation of theinjection above); (3) predicted injection calculation; (4) Minimizationof the instantaneous error; and (5) target tracking is not available.

Open-loop control: patient operated pump: Open-loop control refers tomoderately frequent potentially irregular actuation of insulin treatmentvia an insulin pump capable of delivering basal rate and boluses, wherethe regulator does not control the pump automatically but insteadadvises the patient on the amount and the timing of insulin doses. Forirregular-actuation implementation a monitoring system should be addedto the observer to alarm the user that action is needed, for example ifthe predicted glucose concentration does not come back to target withinreasonable time. (1) observer triggered by arriving sensor data and notonly by actuation time; (2) modified injection for boluses (seecalculation of the injection above); (3) subtraction of basal rate frominjection scheme; (4) possibility of regular actuation schedule ofseveral hours (during the day); (5) predicted injection calculation; and(6) minimization of the instantaneous error.

Closed-loop control refers to insulin treatment via an insulin pumpwhere actuation is automatically transmitted from the algorithm to thepump. It this particular case actuation is regular and frequent. Theuser maintains oversight over the insulin amount injected. In thisparticular application all the proposed modules are applicable, and theactuation frequency can be optimized for minimal energy expenditure ofthe pump. Simulation provides insight into such control as shown in FIG.5, which shows a simulation of automatic control using SC sensor with 1minute actuation time, standard Kalman Filter, predictive actuation.

Although the present invention has been described in considerable detailwith reference to certain preferred versions thereof, other versions arepossible. Therefore, the spirit and scope of the appended claims shouldnot be limited to the description of the preferred versions containedherein.

The reader's attention is directed to all papers and documents which arefiled concurrently with this specification and which are open to publicinspection with this specification, and the contents of all such papersand documents are incorporated herein by reference.

All the features disclosed in this specification (including anyaccompanying claims, abstract, and drawings) may be replaced byalternative features serving the same, equivalent or similar purpose,unless expressly stated otherwise. Thus, unless expressly statedotherwise, each feature disclosed is one example only of a genericseries of equivalent or similar features.

Any element in a claim that does not explicitly state “means for”performing a specified function, or “step for” performing a specificfunction, is not to be interpreted as a “means” or “step” clause asspecified in 35 U.S.C § 112, sixth paragraph. In particular, the use of“step of” in the claims herein is not intended to invoke the provisionsof 35 U.S.C § 112, sixth paragraph.

What is claimed is:
 1. A method for determining an insulin dosingrecommendation comprising: applying an adaptive observer scheme toestimate a present physiological state of a patient based on at leastone measurement of blood glucose concentration in the patient, whereinthe adaptive observer scheme models the present physiological state ofthe patient using a dynamic model; and applying the estimatedphysiological state to a Linear Quadratic regulator to determine aninsulin dosing recommendation based on the present physiological stateso estimated.